Description des positions polymorphes

library(dplyr)
library(ggplot2)
library(viridis)
library(readr)
library(magrittr)

col_adn <- function() col_factor(c("a", "t", "c", "g", "n"))
col_ADN <- function() col_factor(c("A", "T", "C", "G", "N"))

# 10'' environ
data <- readr::read_csv(
  "huge-strong.csv",
  col_types = cols(
    baseCall     = col_adn(),
    calledBase   = col_ADN(),
    uncalledBase = col_ADN(),
    ## areaUncalledBase = col_double(na = c("-1")),
    mutant       = col_factor(c("strong", "weak")),
    A.base       = col_ADN(),
    B.base       = col_ADN(),
    first.call   = col_ADN(),
    second.call  = col_ADN(),
    second.area  = col_double(),
    snp_         = col_factor(c("snp"))),
  na = c("-", "", "-1")) %>%
  select(
    plas = plasmide, # nom de la construction
    mutt = mutant, # weak ou strong
    base = baseCall, # 
    ubas = uncalledBase, # base non appelée
    fcal = first.call,
    scal = second.call,
    qpos = SCF.position, # position sur la séquence
    recv = A.base,
    donn = B.base,
    air1 = area, # aire du pic primaire
    air2 = areaUncalledBase,
    rel1 = relativeAreaCalledPeak, # 
    rel2 = relativeAreaUncalledBase,
    refp = reference.position, 
    ## upos = positionUncalledBase, # position de la base non appelée sur le spectrogramme 
    ## rel2 = parse_double(relativeAreaUncalledBase, na = c("-1", NA)),
    ## file = file,
    ## anls = analysis,
    ## datm = datum,
    ## fara = first.area,
    ## sara = second.area,
    rat1 = uncorrected.proportion.A,
    ## rat2 = uncorrected.proportion.B,
    snp  = snp_,
    qual = quality,
    comt = comments,
    spos = position, # position sur le spectrogramme
    ## spos1= sposition, # position sur le spectrogramme
    A    = A,
    C    = C,
    G    = G,
    T    = T
  ) %>%
  mutate(base = toupper(base),
         plas = gsub("-1073.ab1", "", plas),
         rel2 = parse_double(.$rel2, na = c("-1")),
         air2 = parse_double(.$air2, na = c("-1")))

## réaffecte les niveaux de facteur. 
##
## | facteur | commentaire                                           |
## |---------|-------------------------------------------------------|
## | D       | WARNING: bases could not be matched                   |
## | C       | WARNING: multiple peak span, data are not trustworthy |
## | B       | WARNING: primary peak does not match the reference    |
## | A       | processed normally                                    |

## il n'y a pas dans les données de facteur C. d'où les trois niveaux de
## facteur. sensible aux variations de paramètre polySNP à mon avis. à
## surveiller.
data$comt <- factor(data$comt)
data$comt <- plyr::mapvalues(data$comt,
                             from = levels(data$comt),
                             to = c("D", "B", "A"))

##                                        # algo :
## par plasmide,
## trouve la première position alignée avec la référence
## ajuste les positions avant celle là en fonction
## ajuste les positions après en fonction. 
##
## stratégie : faire une table qui associe par position la position de référence
## connue la plus proche, inner_join on it, calcul de la différence entre la
## position de référence la plus proche et la position actuelle

data %>%
  filter(!is.na(refp)) %>% # ne garde que les positions ou la ref est connue
  mutate(diff = refp - qpos) %>% # calcule la difference entre la position sur
                                 # la sequnece et la position sur la ref
  select(plas, diff) %>% # ne garde que les colonnes avec le nom du plasmide et
                         # la différence
  group_by(plas) %>%  # par plasmide
  summarise(diff = diff %>% mean() %>% round() ) %>% # compte la moyenne des
                                                     # différences, arrondie à
                                                     # l'entier le plus proche
  inner_join(data , ., by = "plas") %>% # regroupe avec la table mère
  mutate(apos = qpos + diff) -> # calcule le décalage entre la position sur la
                                # séquence et la position sur la référnece
  data


theme_set(
  theme_minimal(base_size = 10, base_family = "Courier") +
  theme(panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.text = element_text(face = "bold", size = 12)))


## rajoute une colonne longueur qui détermine la longueur de la trace de
## conversion par plasmide.
data %>%
  rowwise() %>%
  filter(fcal == donn) %>%
  ungroup() %>%
  group_by(plas, mutt) %>%
  summarise(len = max(refp) - min(refp)) %>%
  inner_join(x = data, y = .) ->
  data



# #' Cette fonction retourne un graphique ggplot2, avec en ordonnée les
# #' constructions et en abscisse la position sur la séquence de référence.
# #' Chaque point représente un pic ou phred a détecté un pic secondaire. Les
# #' axes verticaux représentent les positions de snp.
# #' @title plot_second_peak
# #' @param data la table de donnée
# #' @param mutt "strong" ou "weak", ne représente que les positions pour les
# #'   constructions weak ou strong
# #' @param length classe les séquences en fonction de leur longueur. 
# #' @param snp_only ne représente que les positions de SNP introduits. 
# #' @param clean enlève les séquences qui ont trop de polymorphes. 
# #' @return un graphique ggplot
plot_second_peak <-
  function(data, sw = NULL, length = FALSE, clean = TRUE, snp_only = FALSE)
{
  stopifnot(
    sw == "weak" || sw == "strong" || is.null(sw),
    typeof(length)   == "logical",
    typeof(snp_only) == "logical" 
  )

  find_length <-
    function(data) data %>% arrange(len) %>% {.$plas %>% unique()}

  find_snppos <-
    function(data) data %>%
                     filter(!is.na(snp)) %>%
                     { if (!is.null(sw)) filter(., mutt == sw) else . } %>%
                     group_by(refp) %>%
                     summarise(count = n()) %>% {.$refp}

  data %>%
    # si clean is true, ne garde que les séquences absentes de seqcrade
    { if (clean) filter(., !(plas %in% seqcrade)) else . } %>%  
    # filter pour ne garder que les séquneces strong ou weak
    { if (!is.null(sw))   filter(., mutt == sw) else . } %>%
    # filtrer pour ne garder que les snp
    { if (snp_only)       filter(., !is.na(snp))  else . } %>%
    filter(ubas != "N") %>%
    ggplot(
      data = .,
      aes(x = apos,
          # si length, trier par longueur de conversion tract. 
          y = { if (length) (factor(plas, levels = find_length(data)))
                else plas },
          color = ubas, fill = ubas, label = ubas)
    ) +
    geom_vline(xintercept = find_snppos(data), color = "gray", alpha = 0.5) +
    geom_point(alpha = 0.6) +
    xlab("Position sur la sequence de reference") +
    ylab("Transformant") +
    theme(legend.position = "bottom") +
    ## scale_color_(discrete = TRUE, name = "Base\nSecondaire") +
    scale_color_brewer(name = "Base\nSecondaire", palette = "Set1") +
    scale_fill_discrete( name = "Base\nSecondaire")
}

## #' cette fonction définit si une base est strong ou weak. Si la base en question
## #' est NA, alors renvoit un charactere vide.
## #' @param base : une base d'ADN.
base_caller <- function(base)
{
  if (is.na(base) || base == "N" ) { "" }
  else if (base == "A" || base == "T") { "W" }
  else { "S" }
}


find_length <-
  function(data) data %>% arrange(len) %>% {.$plas %>% unique()}

plot_first_second <-  function(data, sw = NULL, snp_only = FALSE) {

  find_snppos <-
    function(data) data %>%
                     filter(!is.na(snp)) %>%
                     { if (!is.null(sw)) filter(., mutt == sw) else . } %>%
                     group_by(refp) %>%
                     summarise(count = n()) %>% {.$refp}

  data %>%
    filter(ubas != "N") %>%
    filter(!(plas %in% seqcrade)) %>%
    { if(!is.null(sw)) filter(., mutt == sw) else . } %>%
    { if(snp_only) filter(., snp == "snp")   else . } %>%
    rowwise() %>%
    mutate(type1  = base_caller(base),
           type2  = base_caller(ubas)) %>%
    ungroup() %>%
    ggplot(aes(x     = apos,
               y     = factor(plas, levels = find_length(.)),
               label = paste0(type1, type2),
               fill  = paste0(type1, type2))) +
    ## geom_tile(alpha = 0.2) +
    geom_point(aes(color = factor(paste0(type1, type2))), alpha = 0.5, size = 4) +
    geom_text(aes(
      color = factor(paste0(type1, type2))),
      size = 2,
      ## alpha = 0.6,
      family = "Courier",
      fontface = "bold") +
    geom_vline(xintercept = find_snppos(data),
               color = "gray", alpha = .5
               ) +
    scale_fill_brewer(palette = "Set1") +
    scale_color_brewer(palette = "Set1") +
    xlab("Position sur la sequence de reference") +
  ylab("") +
  ggtitle(paste("Constructions", toupper(sw))) +
  guides(colour = guide_legend(title = ""), fill = "none") +
    theme(legend.position = "bottom",
          panel.grid.major.y = element_line(color = "gray", linetype = "dotted"))
}

plot_snp_pos <- function(data, sw) {

  faclevel <-
    function(x, y) factor(paste0(x, y), levels = c("SS", "SW", "WS", "WW", "S", "W"))

  data %>%
    filter(snp == "snp") %>%
    filter(mutt == sw) %>%
    filter(!(plas %in% seqcrade)) %>%
    rowwise() %>%
    mutate(type1 = base_caller(base),
           type2 = base_caller(ubas)) %>%
    ungroup() %>%
    ggplot(aes(x = apos,
               y = factor(plas, levels = find_length(.)),
               label = faclevel(type1, type2),
               ## fill = faclevel(type1, type2),
               color = faclevel(type1, type2))) +
    geom_point(alpha = 0.4, size = 4) +
    geom_text(size = 2, family = "Courier", fontface = "bold") +
    ## scale_fill_brewer(palette = "Set1") +
    scale_color_viridis(discrete = TRUE) +
    labs(x = "Position sur la sequence de reference", y = "") +
    ggtitle(paste("Transformants", toupper(sw))) +
    guides(colour = guide_legend(title = ""), fill = "none") +
    theme(legend.position = "bottom",
          panel.grid.major.y = element_line(color = "gray", linetype = "dotted"))

}

Analyses des positions polymorphes

Toutes les positions

Le graphique suivant représente toutes les positions par séquence auxquelles phred détecte un pic secondaire. La couleur représente la base au pic secondaire. Les traits verticaux représentent les positions de SNP introduites. Le trait horizontal démarque les donneurs STRONG (en bas) des WEAK.

plot_second_peak(data, clean = FALSE) +
  geom_hline(yintercept = 84.5, linetype = "dotted")

Certaines séquences sont extrêmement bruitées. Le graphique suivant représente le nombre de positions polymorphe par plasmide.

data %>%
  filter(ubas != "N") %>% # ne regarde que les positions où on a du
                          # polymorphisme
  group_by(plas,mutt) %>%
  summarise(ucnt = n()) %T>% # par plasmide, compte le nombre de positions
                             # polymorphes
  { print( ggplot(data = .) +
           geom_point(aes(ucnt, plas, color = mutt)) +
           geom_vline(xintercept = 80, color = "red", linetype = "dotted") +
           scale_color_brewer(palette = "Set1", guide = FALSE)
          )} %>%
  filter(ucnt > 80) %>% # au vu de la distribution, 80 semble un bon cutoff de
                         # propreté
  {.$plas} ->
  seqcrade # sortie dans un vecteur.

J’ai utilisé le seuil de 80 positions polymorphes pour discriminer les séquences propres des séquences “sales” (cutoff en rouge). Sans ces même séquences :

plot_second_peak(data, clean = TRUE) +
  geom_hline(yintercept = 79.5, linetype = "dotted")

## #' Seulement les séquences des transformants _strong_
## #+ seqstrong, fig.height = 7.5
## plot_second_peak(data, sw = "strong")

## #' Seulement les séquences des transformants _weak_
## #+ seqweak, fig.height = 7.5
## plot_second_peak(data, sw = "weak")

Si on ne regarde que les positions de SNP

plot_second_peak(data, snp_only = TRUE) + geom_hline(yintercept = 11.5, linetype = "dotted")

On voit plusieurs choses :

1. il y a nettement plus de positions polymorphes dans nos marqueurs avec les donneurs WEAK qu’avec les donneurs STRONG (184 contre 26),

data %>%
  filter(!is.na(snp)) %>%
  filter(!(plas %in% seqcrade)) %>%
  filter(ubas != "N") %>%
  group_by(mutt) %>%
  summarise(count = n()) %>%
  ggplot(data = ., aes(x = mutt, y = count, fill = mutt )) +
  geom_bar(stat = "identity", width = 0.01) +
  geom_point(aes(color = mutt)) +
  scale_fill_brewer(palette = "Set1" , guide = FALSE) +
  scale_color_brewer(palette = "Set1", guide = FALSE) +
  coord_flip() +
  labs(x = "Nombre de\npositions polymorphes", y = "Donneurs")

2. le nombre de positions polymorphes par séquence est assez variable.

data %>%
  filter(!is.na(snp)) %>%
  filter(!(plas %in% seqcrade)) %>%
  filter(ubas != "N") %>%
  group_by(mutt, plas) %>%
  summarise(count = n()) %>%
  ungroup() %>%
  arrange(desc(count)) %>%
  ggplot(data = .,
         aes(y = count,
             x = factor(plas, levels = {.$plas}),
             fill = mutt, color = mutt)) +
  geom_point() +
  geom_bar(stat = "identity", width = 0.1, alpha = 0.3) +
  coord_flip() +
  scale_fill_brewer(palette = "Set1", guide = FALSE) +
  scale_color_brewer(palette = "Set1", guide = FALSE) +
  labs(x = "Transformants", y = "Nombre de positions polymorphes")

3. les pics secondaires ne sont pas nécessairement spécifiques de nos marqueurs :

data %>%
  filter(ubas != "N") %>%
  filter(!(plas %in% seqcrade)) %>%
  group_by(snp) %>%
  summarise(count = n()) %>%
  ggplot(data = .,
         aes(x = snp, y = count)) +
  geom_bar(stat = "identity", width = 0.05, fill = "gray") +
  geom_point() +
  coord_flip() 

4. mais la fréquence de pic secondaire est plus élevée dans les sites marqueurs que dans les sites normaux.

data %>%
  filter(!(plas %in% seqcrade)) %>%
  filter(ubas != "N") %>%
  group_by(plas) %>%
  summarise(totalpoly = n()) %>%
  inner_join(x = data, y = . ) %>%
  filter(ubas != "N") %>%
  filter(!(plas %in% seqcrade)) %>%
  group_by(mutt, plas, snp, totalpoly) %>%
  summarise(count = n()) %>%
  ungroup() %>%
  mutate(freq = count / totalpoly) %>%
  filter(!is.na(snp)) %>%
  arrange(desc(freq)) %>%
  ggplot(data = .,
         aes(x = factor(plas, levels = {.$plas}),
             y = freq)) +
  geom_bar(aes(fill = mutt), stat = "identity", width = 0.2) +
  geom_point(aes(size = count, color = count)) +
  coord_flip() +
  scale_color_viridis(begin = 0.2, end = 0.9,
                      name = "Nombre de positions\navec pic secondaire",
                      guide = guide_legend(direction = "horizontal",
                                           label.position = "bottom")) +
  scale_fill_viridis(discrete = TRUE, begin = 0.8, end = 0.2,
                     name = "Type de donneur",
                     guide = guide_legend(direction = "horizontal")) +
  scale_size_area(max_size = 10, guide= FALSE) +
  labs(x = "Donneur", y = "Frequence des sites avec\npic secondaire correspondant aux marqueurs") +
  theme(legend.position = c(0.8, 0.8))

data %>%
  filter(ubas != "N") %>%
  group_by(plas) %>%
  summarise(totalpoly = n()) %>%
  inner_join(x = data, y = . ) %>%
  filter(ubas != "N") %>%
  filter(!(plas %in% seqcrade)) %>%
  group_by(mutt, plas, snp, totalpoly) %>%
  summarise(count = n()) %>%
  ungroup() %>%
  mutate(freq = count / totalpoly) %>%
  arrange(desc(totalpoly)) %>%
  mutate(snp = factor(ifelse(is.na(snp), "autre", "snp"),
                      levels = c("autre", "snp"))) %>%
  ggplot(data = ., aes(x = factor(plas, levels= {.$plas}), y = count )) +
  geom_bar(aes(fill = snp), stat = "identity", width = 0.3, position = "dodge") +
  geom_point(aes(color = snp)) +
  scale_color_viridis(begin = 0.2, end = 0.8, discrete = TRUE ) +
  scale_fill_viridis(begin = 0.2, end = 0.8, discrete = TRUE ) +
  geom_hline(yintercept = 0, color = "gray") +
  labs(x = "Donneur", y = "Nombre de positions") +
  theme(
    axis.text.x = element_text(angle = 90, size = 6),
    panel.grid.major.x = element_line(linetype = "dotted", color = "gray")
  )

J’ai ensuite voulu regarder la répartition par séquence de ces seconds pics. La base majoritaire correspond-elle plus souvent à la base donneuse, à la base receveuse, ou aucun des deux ?

Le graphique suivant représente donc toutes les positions où phred détecte un pic secondaire, encore une fois. La couleur du point représente les différents cas de figure : rouge pour les cooccurences SS, bleu pour les cooccurences SW (pic majoritaire S et pic secondaire W), vert pour l’inverse et violet pour les cas WW.

plot_first_second(data)

Si on compte le nombre de cooccurences différentes :

data %>%
  filter(ubas != "N") %>%
  filter(!(plas %in% seqcrade)) %>%
  rowwise() %>%
  mutate(altern = paste0(base_caller(base), base_caller(ubas)) %>%
           factor(levels = c("SS", "SW", "WS", "WW", "S", "W"))) %>%
  group_by(mutt, altern) %>%
  summarise(count = n()) %>%
  ggplot(data = .,
         aes(x = altern, y = count, fill = altern)) +
  geom_point(aes(color = altern, size = count)) +
  geom_bar(stat = "identity", width = 0.1, alpha = 0.4) +
  facet_grid(mutt~.) +
  scale_size(guide = FALSE) +
  scale_fill_brewer( palette = "Set1", guide = FALSE) +
  scale_color_brewer(palette = "Set1", guide = FALSE) +
  labs(y = "Nombre de positions polymorphes",
       x = "Type de coocurrence") +
  theme(legend.position = "bottom") +
  coord_flip()

Il y a autant de coocurrences WS (donneur / sauvage) quand le donneur est weak que de cooccurrences SW (donneur / sauvage) quand le donneur est strong. Par contre, il y a plus de cooccurrences SW (sauvage / donneur) quand le donneur est weak que de cooccurences WS (sauvage / donneur) quand le donneur est strong. Il y a également un nombre relativement élevé de cooccurences SS et WW.

## #+ firstsecondweak, fig.height = 7
## plot_first_second(data, "weak")

Si on ne regarde que les positions de SNP pour les donneurs WEAK.

plot_first_second(data, "weak", TRUE)

Idem, seulement les SNP pour les donneurs STRONG.

plot_first_second(data, "strong", TRUE)

Certains reads ne sont pas homogènes : tantôt la base donneuse est majoritaire, tantôt la base receveuse. Le graphique suivant est un alignement de toutes les positions de SNP pour les donneurs WEAK, classés par longueur de trace de conversion.

plot_snp_pos(data, "weak") 

Idem pour les transformants strong.

plot_snp_pos(data, "strong")

On en a conclu plusieurs choses :

1. il semble que les couleurs associées soit plutôt bleu avec jaune et violet avec vert. Autrement dit les pics majoritaires dans la trace de conversion sont plutôt weak chez les donneurs weak, et les pics majoritaires en dehors de la trace correspondent à la base strong.

2. Dans l’hypothèse où ce seraient des contaminations : on dispose du plan de plaque, on peut donc essayer de tester l’association entre les contaminations qu’on observe et les séquences présentes dans les transformants des autres puits. Une sorte de démineur des conta. La méthode n’est encore pas établie.

3. On retrouve toujours les bases strong (en vert) au milieu des séquences de weak (en jaune). Ça pourrait être dû à des traces de conversions complexes. En tout les cas ces mutations ne sont pas associées à des seconds pics, pour 4 des 6 cas que Laurent avait détecté auparavant.

4. Concernant le problème des contaminations : on s’est dit qu’il pouvait subsister l’ADN du plasmide, adsorbé sur la bactérie transformante et qu’on isole sur le milieu. La bactérie isolée, en cours de croissance, pourrait utiliser cet ADN pour un deuxième évènement de recombinaison. Si au lieu d’un évènement de recombinaison, on a plusieurs évènements, qui utilisent cet ADN résiduel comme matrice, ça pourrait expliquer les profils de polymorphismes qu’on observe, qui semblent associé (globalement) spécifiquement à nos SNP. On veut donc traiter les clones, avant d’étaler, par de la DNAse, pour éliminer les traces de plasmides résiduels. # distribution de la longueur de trace de conversion

data %>%
  filter(!is.na(snp)) %>% # ne regarde que les positinos de snp
  filter(!(plas %in% seqcrade)) %>% # enleve les sequences de faible qualité
  rowwise %>% #
  filter( base == donn ) %>% # détermine la première base sur la référence qui
                             # n'est pas une base donneuse.
  ungroup() %>%
  group_by(plas, mutt) %>%
  summarise(switchp = min(refp, na.rm = TRUE) ) %>% # appelée switchp
  group_by(mutt, switchp) %>%
  summarise(count = n()) %>% # détermine la distribution par mutant des ces
                             # switchp
  ggplot(., aes(switchp, count)) +
  geom_bar(aes(fill = mutt), stat = "identity",
           alpha = 0.5, position = "dodge", width = 0.5) +
  geom_point(aes(color = mutt)) +
  scale_fill_discrete(guide = FALSE) +
  scale_color_brewer(palette = "Set1", name = "Donneur") +
  labs(x = "Position de basculement")

On ne voit pas une variation affolante entre la longueur des traces de conversions issues des donneurs weak et strong. # Variations des scores de confiance

Chaque fois qu’un pic secondaire est détecté, un ratio d’aire est calculé entre le pic primaire et le pic secondaire (ratio = 1 => pics équivalents). J’ai voulu voir la façon dont ces ratios étaient distribués sur les séquences. On peut regarder cette variabilité intra et inter-séquence.

Variabilité intra séquence

Le graphique suivant représente la distribution des scores de confiance par read. La couleur et la taille des points correspondent au ratio des pics (secondaire / primaire, toujours < 1).

data %>%
  filter(ubas != "N") %>%
  ## filter(mutt == "weak") %>%
  filter(!(plas %in% seqcrade)) %>%
  ggplot(data = .,
         aes(x = apos, y = plas)) +
  geom_point(aes(size = rel2, color = rel2), alpha = 0.8) +
  scale_color_viridis(begin = 1, end = 0, name = "Ratio") +
  scale_size_area(max_size = 3, guide = FALSE) +
  labs(x = "Position sur la reference", y = "Transformants") +
  theme(legend.position = "bottom",
        panel.grid.major.y = element_line(color = "gray", linetype = "dotted"))

Il semble que les scores varient assez peu au sein d’une même séquence. C’est ce que le graphique suivant cherche à montrer.

data %>%
  filter(ubas != "N") %>%
  ## filter(mutt == "strong") %>%
  filter(!(plas %in% seqcrade)) %>%
  ggplot(data = ., aes(x = apos, y = rel2)) +
  geom_point(aes(size = rel2, color = rel2), alpha = 0.8) +
  geom_line(aes(group = plas, color = rel2), alpha = 0.3) +
  geom_vline(xintercept = data %>% filter(!is.na(snp)) %>% group_by(refp) %>% summarise(count = n()) %>% {.$refp},
             color = "gray", size = 0.1) +
  scale_color_viridis(begin = 1, end = 0, guide = FALSE) +
  scale_size_area(max_size = 3, guide = FALSE) +
  labs(x = "Position sur la reference", y = "Ratio\nPic Secondaire /\nPic Primaire") +
  theme(legend.position = "bottom",
        panel.grid.major.y = element_line(color = "gray", linetype = "dotted"))

Hormis quelques points, les lignes se croisent assez peu, d’une position à l’autre. Le graphique suivant représente la mediane et la déviation absolue à la médiane des ratios de pic, globalement assez ressérés autour de 0.15. 0.2 est un score assez conservateur, il correspond clairement à un pic secondaire, lorsqu’on le vérifie visuellement sur les spectrogrammes.

data %>%
  filter(ubas != "N") %>%
  filter(!(plas %in% seqcrade)) %>%
  group_by(plas, mutt) %>%
  summarise(med = median(rel2, na.rm = TRUE),
            mad = mad(rel2, na.rm = TRUE)) %>%
  arrange(mad) %>%
  ggplot(data = .,
         aes(y = plas, x = med, color = med )) +
  geom_point() +
  geom_errorbarh(aes(xmax = med + mad, xmin = med - mad )) +
  scale_color_viridis(begin = 1, end = 0, guide = FALSE) +
  theme(legend.position = "bottom") +
  labs(x = "Médiane des ratios par read")

J’ai voulu regarder sur l’ensemble des séquences, la fréquence des sites avec un pic secondaire par position. L’idée est de montrer que les pics secondaires sont plus fréquents aux positions de SNP maqueurs qu’on a introduit qu’aux autres positions.

data %>%
  mutate(has_sp = ifelse(ubas != "N", TRUE, FALSE)) %>%
  filter(has_sp) %>%
  group_by(apos) %>%
  summarise(n_ps = n()) %>%
  inner_join(x = data, y = .) %>%
  group_by(apos, n_ps) %>%
  summarise(n_apos = n()) %>%
  mutate(freq = n_ps / n_apos) %>%
  ggplot(data = ., aes(apos, freq)) +
  geom_point(aes(color = freq, size = freq)) +
  geom_line(alpha = 0.1) +
  scale_x_continuous(breaks =
                       data %>% filter(!is.na(snp), mutt == "weak") %>%
                       group_by(refp) %>%
                       summarise(count = n()) %>% {.$refp} ) +
  scale_y_continuous(limits = c(0, 0.30)) + 
  scale_color_viridis(begin = 0.95, end = 0) +
  ## scale_y_log10() +
  guides(color = "none", alpha = "none", size = "none") +
  ggthemes::theme_fivethirtyeight() +
  xlab("Position sur la reference") +
  ylab("Frequence de pic secondaire") +
  ggtitle("Frequence des bases avec pic\nsecondaire a une position donnee") +
  theme(
    panel.grid.major.y = element_line(linetype = "dotted", color = "gray", size = 0.1),
    panel.grid.major.x = element_line(linetype = "dotted", color = "gray", size = 0.4),
    axis.text.x = element_text(size = 6),
    axis.text.y = element_text(size = 6),
    axis.title.x = element_text(color = "black", size = 8)) 

Type de transitions

J’ai enfin voulu regarder quelles types de transitions étaient favorisées selon les donneurs.

purpyr_caller <- function(base)
{
  if (is.na(base) || base == "N" ) { "N" }
  else if (base == "A" || base == "G") { "pur" }
  else { "pyr" }
}

## data %>%
##   filter(!is.na(snp)) %>%
##   filter(base != recv) %>%
##   mutate(trans = paste0(base, recv)) %>%
##   group_by(mutt, trans) %>%
##   summarise(count = n()) %>%
##   ggplot(data = .,
##          aes(x = trans, y = count )) +
##   geom_point() +
##   geom_bar(aes(fill = mutt), stat = "identity")

data %>%
  filter(!is.na(snp)) %>%
  filter(base != recv) %>%
  rowwise() %>%
  mutate(purpyr_base = purpyr_caller(base),
         purpyr_wt   = purpyr_caller(recv),
         purpyr_trans= paste0(purpyr_wt,"->",purpyr_base)) %>%
  filter(purpyr_base != "N") %>%
  ungroup() %>%
  group_by(mutt, purpyr_trans) %>%
  summarise(count = n()) %>%
  ggplot(data = .,
         aes(x = purpyr_trans, y = count, size = count, color = mutt )) +
  geom_point() +
  scale_color_viridis(begin = 0.1, end = 0.9, discrete = TRUE) +
  scale_size(guide = FALSE) +
  labs(x = "", y = "Nombre de transitions") +
  coord_flip() +
  theme(panel.grid.major.y = element_line(linetype = "dotted", color = "gray"))

Analyse de la répartition des séquences avec pics secondaires dans la plaque

Le but est de faire un plan de plaque pour étudier l’association entre le nombre de position polymorphe et la position dans la plaque.

plaque <- matrix(1:96, nrow = 8, ncol = 12) 
colnames(plaque) <- LETTERS[1:12] 
rownames(plaque) <- letters[1:8]  
plaque <- reshape2::melt(plaque, value.name = "puit")
colnames(plaque) <- c("colonne", "ligne", "puit")

data %>%
  filter(!(plas %in% seqcrade )) %>% # ne garde que les séquences propres 
  filter(!is.na(snp)) %>% # ne garde que les positions de SNP
  filter(ubas != "N") %>% # ne garde que les postions de pic secondaire
  rowwise() %>% 
  mutate(plas = gsub(pat= "p.", rep = "", x = plas) ) %>% # remplace les
                                                          # caractères inutiles
  ungroup() %>%
  group_by(plas, mutt) %>% # compte le nombre d'occurence par transformant
  summarise(count = n()) %>%
  ungroup() %>%
  mutate(plas  = as.integer(plas)) %>% # même types pour le join
  left_join(plaque, y = ., by = c("puit" = "plas")) %>% # joint sur le puit
  mutate(count = replace(count, is.na(count), 0))  %>% # remplace les NA par 0
  filter(!is.na(mutt)) %>% # enlève les positions où on n'a rien
  ggplot(aes(y = factor(colonne, levels = c("h", "g", "f", "e", "d", "c", "b", "a")),
             x = ligne)) +
  geom_point(aes(color = count, size = count)) +
  facet_grid(mutt~.) +
  scale_color_viridis() +
  scale_size_area(max_size = 20) +
  guides(
    colour = guide_legend(title = "Nombre de positions\navec pic secondaire",
                          size = 2),
    size = "none") +
  ggthemes::theme_fivethirtyeight() +
  theme(
    panel.grid.major.y = element_line(linetype = "dotted", color = "gray"),
    panel.grid.major.x = element_line(linetype = "dotted", color = "gray"),
    legend.title = element_text(size = 8),
    plot.title = element_text(size = 18),
    panel.margin = unit(1, "cm"),
    strip.text = element_text(size = 14, face = "bold")
  ) +
  ggtitle("Emplacement des transformants\nmontrant des pics secondaires")

Il semble quand on regarde vaguement comme ça les séquences “contaminées”, ie avec des pics secondaires, qu’elles soient relativement clusterisées autour de quelques positions. On a donc voulu tester par une simulation si la répartition est aléatoire ou non. Si elle ne l’est pas, ça peut vouloir indiquer une contamination.

# reconstruit la plaque weak en matrice
data %>%
  filter(!(plas %in% seqcrade )) %>% # ne garde que les séquences propres 
  filter(!is.na(snp)) %>% # ne garde que les positions de SNP
  filter(mutt == "weak") %>%
  filter(ubas != "N") %>% # ne garde que les postions de pic secondaire
  rowwise() %>% 
  mutate(plas = gsub(pat= "p.", rep = "", x = plas) ) %>% # remplace les
                                                          # caractères inutiles
  ungroup() %>%
  group_by(plas) %>% # compte le nombre d'occurence par transformant
  summarise(count = n()) %>%
  ungroup() %>%
  mutate(plas  = as.integer(plas)) %>% # même types pour le join
  left_join(plaque, y = ., by = c("puit" = "plas")) %>% # joint sur le puit
  mutate(count = replace(count, is.na(count), 0))  %>% # remplace les NA par 0
  {.$count %>% matrix(nrow = 8, ncol = 12)} ->
  weak_plaque

print(weak_plaque)

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,]    0    1    0    0    0    0    2    0    0     0     0     1
## [2,]    0    0    0    0    0    0    0    0    0     7     0     0
## [3,]    0    0    2    0   17    0    0    0    0     1    15     0
## [4,]    3    0    3    0    4    4    0    1   15     2     1     0
## [5,]   15    1    6    0    1    0    0    0   21     2     9     0
## [6,]    0    0    7   10    8    0    0    0    5     0     0     6
## [7,]    0    1    1    1    3    1    0    0    0     0     0     0
## [8,]    0    0    1    0    0    6    0    0    0     0     0     0

# reconstruit la plaque strong en matrice
data %>%
  filter(!(plas %in% seqcrade )) %>% # ne garde que les séquences propres 
  filter(!is.na(snp)) %>% # ne garde que les positions de SNP
  filter(mutt == "strong") %>%
  filter(ubas != "N") %>% # ne garde que les postions de pic secondaire
  rowwise() %>% 
  mutate(plas = gsub(pat= "p.", rep = "", x = plas) ) %>% # remplace les
                                                          # caractères inutiles
  ungroup() %>%
  group_by(plas) %>% # compte le nombre d'occurence par transformant
  summarise(count = n()) %>%
  ungroup() %>%
  mutate(plas  = as.integer(plas)) %>% # même types pour le join
  left_join(plaque, y = ., by = c("puit" = "plas")) %>% # joint sur le puit
  mutate(count = replace(count, is.na(count), 0))  %>% # remplace les NA par 0
  {.$count %>% matrix(nrow = 8, ncol = 12)} ->
  strong_plaque

print(strong_plaque)

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,]    0    0    0    0    0    0    0    0    0     1     0     0
## [2,]    0    0    1    0    0    0    0    0    0     0     0     0
## [3,]    0    0    0    0    0    2    0    0    0     0     0     0
## [4,]    0    0    0    0    0    0    0    0    1     0     0     0
## [5,]    0    0    0    0    1    0    0    0    0     0     0     0
## [6,]    0    0    0    0    0    0    0    0    0     0     0     0
## [7,]    0    0    1    0    0    0    0    0    1     0     0     0
## [8,]    0    1    0    0    1    0    0    0    1     0     0    15

Il semblerait que les séquences contaminées soient réparties uniformément dans les séquences strong, mais qu’elles soient clusterisées autour de deux “secteur” dans la plaque weak. On a vu avec Laurent pour mettre au point une simulation. Celle ci nous permettrait de voir la répartition des puits avec voisins contaminés sous l’hypothèse d’une répartition aléatoire. Si la moyenne du nombre de puits avec voisins contaminés n’est obtenue qu’un petit nombre de fois dans la simulation, alors au risque égal à \(10^{-X}\), avec X le nombre de répétitions, on peut dire que la répartition des puits contaminés n’est pas dûe au hasard.

library(assertthat) # tests

# une fonction qui génère un sampler entre 1 et n.
random_well <- function(n) sample(1:n, 1)

# une fonction qui initialise la plaque à 0 conta par puit
initiate_plate <- function() matrix(rep(0, 96), nrow = 8, ncol = 12)

# une fonction qui transforme l'entrée en matrice. par défaut, crée une plaque
# vide. doublon. conservée par inertie. et flemme.
as.plate <- function(mat = initiate_plate()) matrix(mat, nrow = 8, ncol = 12)

# génère une plaque avec n puits contaminés dedans.
plate_generator <- function(plate = as.plate(), n) {
  assert_that(is.matrix(plate), is.numeric(n), n < 96)

  # une fonction qui contamine un puit dans la plaque s'il ne l'est pas déjà.
  contaminer <- function(mat) {
    assert_that(is.matrix(mat))

    x <- random_well(8)  # tire une ligne au hasard
    y <- random_well(12) # tire une colonne au hasard

    # si le puit n'est pas déjà contaminé, contamine le.
    # sinon, tire un autre puit et contamine le.
    if (mat[x, y] != 1) { mat[x, y] <- 1; return(as.matrix(mat)) }
    else { contaminer(mat = as.matrix(mat)) }
  } 

  # si la plaque n'a pas déjà n puits contaminés, contamine une plaque avec n-1
  # puits, et ainsi de suite.
  # si la plaque a déjà n puits contaminés, renvoit là. 
  if (sum(plate < n )) { plate_generator(as.plate(plate) %>% contaminer(), n - 1) }
  else                 { return(as.plate(plate)) }
}

# génére n plaques contaminées avec n_conta inside.
n_plate_generator <- function(n_plate, n_conta)
  lapply(1:n_plate, function(i) plate_generator(n=n_conta))

### example use :
### ============= 
## n_plate_generator(3, 3)

# une fonction qui teste si les puits voisins sont contaminés. Elle est
# profondément testée. DO NOT EVEN THINK OF TOUCHING IT.
has_nearby_conta <- function(plate, ligne, colonne) {
  assert_that(is.matrix(plate), is.numeric(ligne), is.numeric(colonne))

  x      <- ligne
  y      <- colonne
  p      <- plate
  nearby <- c()

  bord_gauche     <- function(col) ifelse(col == 1 , TRUE, FALSE)
  bord_droite     <- function(col) ifelse(col == 12, TRUE, FALSE)
  bord_sup        <- function(row) ifelse(row == 1 , TRUE, FALSE)
  bord_inf        <- function(row) ifelse(row == 8 , TRUE, FALSE)

  coin_sup_gauche <- function(row, col) ifelse(bord_gauche(col) & bord_sup(row), TRUE, FALSE)
  coin_inf_gauche <- function(row, col) ifelse(bord_gauche(col) & bord_inf(row), TRUE, FALSE)
  coin_inf_droite <- function(row, col) ifelse(bord_droite(col) & bord_inf(row), TRUE, FALSE)
  coin_sup_droite <- function(row, col) ifelse(bord_droite(col) & bord_sup(row), TRUE, FALSE)

  if      (coin_sup_gauche(x, y)) { nearby <- c(nearby, p[x, y+1], p[x+1, y+1], p[x+1, y])}
  else if (coin_inf_gauche(x, y)) { nearby <- c(nearby, p[x-1, y], p[x-1, y+1], p[x, y+1])}
  else if (coin_sup_droite(x, y)) { nearby <- c(nearby, p[x, y-1], p[x+1, y-1], p[x+1, y])}
  else if (coin_inf_droite(x, y)) { nearby <- c(nearby, p[x, y-1], p[x-1, y-1], p[x-1, y])}
  else if (bord_gauche(y))        { nearby <- c(nearby, p[x-1, y], p[x-1, y+1], p[x, y+1], p[x+1, y+1], p[x+1, y])}
  else if (bord_droite(y))        { nearby <- c(nearby, p[x-1, y], p[x-1, y-1], p[x, y-1], p[x+1, y-1], p[x+1, y])}
  else if (bord_sup(x))           { nearby <- c(nearby, p[x, y-1], p[x+1, y-1], p[x+1, y], p[x+1, y+1], p[x, y+1])}
  else if (bord_inf(x))           { nearby <- c(nearby, p[x, y-1], p[x-1, y-1], p[x-1, y], p[x-1, y+1], p[x, y+1])}
  else {
    nearby <- c(nearby,
                p[x-1, y-1], p[x, y-1], p[x+1, y-1],
                p[x-1, y  ],            p[x+1, y  ],
                p[x-1, y+1], p[x, y+1], p[x+1, y+1] )
  }

  ## normalise to 1 if > 1
  nearby[nearby > 0] <- 1
  ## return TRUE if one neighbor is contaminated, false either.
  ## ifelse(sum(nearby) > 0, TRUE, FALSE) 

  ## compte la moyenne du nombre de cas positif pour ce puit 
  return(mean(nearby))
}

# compte la moyenne par plaque de la moyenne du nombre de puits DONT LES
# VOISINS sont contaminés
count_nearby_conta_well <- function(plate) {
  # retrouve les indices des positions par colonne qui sont contaminés
  locate_conta_well <- function(plate) which(plate != 0, arr.ind = TRUE)

  # listes des coordonnées des puits avec conta
  well_list <- locate_conta_well(plate)

  # applique à cette liste la fonction qui permet de retrouver la moyenne du
  # nombre de puit contaminé autour de la position donnée par x[1] et x[2]
  apply(well_list, 1,
        function(x) has_nearby_conta(plate, ligne = x[1], colonne = x[2])) %>%
    # moyenne ces moyennes sur l'ensemble de la plaque
    mean()
}

# compte le nombre de puits contaminés
count_conta <- function(plate) length(which(plate != 0))

count_conta(strong_plaque) # => 11

## [1] 11

count_conta(weak_plaque) # => 35

## [1] 35

# nombre de plaque dans la simulation
n_plaque <- 10000

lance la simulation pour les plaques.

## crée une liste de n_plaque qui contient count_conta(*_plaque) réparties
## aléatoirement sur la plaque.
random_weak        <- n_plate_generator(n_plaque, count_conta(weak_plaque)  ) # takes approx 20''
random_strong      <- n_plate_generator(n_plaque, count_conta(strong_plaque)) # takes approx 20''
random_weak_mean   <- lapply(random_weak  , count_nearby_conta_well) # takes approx 30''
random_strong_mean <- lapply(random_strong, count_nearby_conta_well) # takes approx 30''

# une fonction pour convertir les résultats de lapply en vecteur puis data.frame
# et représenter la distribution.
plot_mean_distrib <- function(rdm_list) {
  assert_that(is.list(rdm_list))

  rdm_list %>% unlist() %>% data.frame(wellcount = .) %>%
    ggplot(aes(wellcount)) + geom_histogram(binwidth = 0.0005, fill = "gray") +
    labs(x = "Moyenne par plaque de la moyenne du\nnombre de puits voisin contamine par puit", y = "") +
    theme(panel.ontop = TRUE,
          panel.grid.major.y = element_line(colour = "white", size = 0.5, linetype = "dotted"),
          panel.grid.minor.y = element_line(colour = "white", size = 0.5, linetype = "dotted"))
}

plot_mean_distrib(random_weak_mean) +
  geom_vline(xintercept = count_nearby_conta_well(weak_plaque), color = "red") 

plot_mean_distrib(random_strong_mean) +
  geom_vline(xintercept = count_nearby_conta_well(strong_plaque), color = "red") 

# en admettant une distribution normale, calcule la statistique associée 
pnorm(q = count_nearby_conta_well(weak_plaque),
      mean = mean(random_weak_mean %>% unlist()),
      sd   =   sd(random_weak_mean %>% unlist()),
      lower.tail = FALSE) %>% print() # => 0.003313558

## [1] 0.003663744

pnorm(q = count_nearby_conta_well(strong_plaque),
      mean = mean(random_strong_mean %>% unlist()),
      sd   =   sd(random_strong_mean %>% unlist()),
      lower.tail = FALSE) %>% print() # => 0.8288228

## [1] 0.8339551